Superoperator nonequilibrium Green's function theory of many-body systems; 
Applications to charge transfer and transport in open junctions 
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Nonequilibrium Green's functions provide a powerful tool for computing the dynamical response 
and particle exchange statistics of coupled quantum systems. We formulate the theory in terms of 
the density matrix in Liouville space and introduce superoperator algebra that greatly simplifies 
the derivation and the physical interpretation of all quantities. Expressions for various observables 
are derived directly in real time in terms of superoperator nonequilibrium Green's functions 
C^l , (SNGF), rather than the artificial time- loop required in Schwinger's Hilbert-space formulation. 

Applications for computing interaction energies, charge densities, average currents, current induced 
fluorescence, electroluminescence and current fluctuation (electron counting) statistics are discussed. 
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I. INTRODUCTION 



C"| \ Nonequilibrium Green's function theory (NEGFT) [H, @, ^ S S B 0] is widely used for computing electron transport 
and optical properties of open many-body systems, such as semiconductors Q, metals Q, molecular wires and scanning 
tunneling microscopy (STM) junctions [icj|. It has also been applied to X-ray photoemission spectroscopy [TTL IT2L [l3l| . 
General fluctuation-dissipation relations may be derived for nonlinear response and fluctuations using this formalism 
[OLE!!- In NEGFT, the diagrammatic perturbative expansions are carried on a closed time-loop which includes two 
branches with forward and backward time propagation, respectively. The closed-time-path was first introduced by 
Schwinger [l| and further developed by Keldysh and others 0, H, [l6| . The physical time t is replaced by a loop-time r 
which runs clockwise on a contour (Fig. la) such that when r goes from tj to t/, the physical time t runs first forward 
and then backward. This technique allows to treat nonequilibrium systems using quantum field theory tools originally 
£h ' developed for equilibrium systems. 

The NEGFT has been remarkably successful in describing a broad range of non-equilibrium systems and phenomena. 
Frequency domain non- linear susceptibilities can be readily interpreted by an evolution on the loop [171 ] . However 
loop approach does not offer an obvious physical real-time intuition for the various quantities and approximations. 
This is only possible at the end of the calculation, when one transforms r to real time. A highly desirable alternative, 
real-time, formalism has been dev elop ed using two time ordering prescriptions known as the "single time" and the 
"physical" representations (l4l. Il5l. \lK fl9| . For reasons that will become clear later, we shall denote these the PTBK 
(Physical-time bra-ket) and the PTSA (Physical-time symmetric-antisymmetric) representations, respectively. We 
shall refer to Schwinger's formulation as the closed time path loop (CTPL). In Ref. [13] NEGFT for bosons was 
formulated in the PTSA representation by considering different evolutions in the forward and backward branches of 
the time-loop. A generating functional for the Keldysh Green's functions was constructed by introducing different 
artificial fields which couple to operators in the forward and the backward branches of the loop. This is similar to the 
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Liouville space [2fJ formulation of NEGF presented in Ref. [2l| where Hedin's equations [22|, [23|, |2J] were generalized 
to an open system. The PTSA representation was finally obtained in Ref. [1J] in two steps; first formulating the 
. ' theory in the PTBK and subsequently making a matrix transformation on the generating functional to obtain the 
^ . PTSA representation. 

In this article we show how by formulating the many-body problem using superoperators in Liouville space, the 
PTSA representation can be used from the outset without introducing any artificial fields. The matrix transformation 
between the two representations is performed on the superoperators themselves rather than on some expectation values 
(Green's functions or generating functions). The time evolution of superoperators in the PTSA representation can then 
be calculated directly from the microscopic equations of motion for superoperatrors, totally avoiding the intermediate 
PTBK representation. The observables are computed directly in terms of the retarded, advanced and correlation 
superoperator nonequilibrium Green's functions (SNGF). For completeness, we introduce bosonic superoperators and 
outline the superoperator approach for bosons. 

Applications are made to inelastic resonances in STM currents [25] and current induced fluorescence in STM 
junctions [261 ]. A key ingredient of the SNGF approach is a simple time-ordering prescription of superoperators that 
provides an intuitive and powerful bookkeeping device for all interactions, and maintains a physical picture based on 
the density-matrix. The physical significance of the various Green's becomes obvious. 

In contrast to Hilbert space CTPL which targets the wave function, the Liouville space PTBK and PTSA formu- 
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lations are based on the density matrix. The backward evolution in Hilbert space is replaced by the simultaneous 
evolution of the bra and ket of the density matrix, which can have a different evolution (but always forward in time!) 
in Liouville space (Fig. lb). Our notation appears naturally when using the density matrix. The "single time", 
PTBK, works with superoperators which act on the bra and the ket. In the "physical", PTSA, representation which 
most directly resembles some observables, we work with linear combinations (sums and differences) of the PTBK 
superoperators. 

The SNGF has been first developed and applied to closed interacting many-bod y sy stems; nonlinear optical response 
of excitons in molecular aggregates [13] and non-equilibrium Van der Waals forces |29l . [30j . Properties of the interacting 
system were determined by the density response as well as the correlated densityfluctuations [29l . l3lj of the individual 
sub-systems. The same approach has been used to resolve the causality paradox [32l l33j of density-functional theory [3^, 
l35j . Here we extend this method to open many-body systems with overlapping charge densities so that electron transfer 
is possible and the number of particles in each subsystem can fluctuate. The SNGF appear naturally as a consequence 
of the time ordering of the ket and the bra evolution of the density matrix which gives rise to Liouville space pathways 
(LSP)[20]. Each SNGF represents a particular combination of LSPs. This establishes a connection between many- 
body nonequilibrium theory and the standard formulation of nonlinear optical response (20l | and paves the way for 
computing nonlinear optical properties of quantum junctions. 

Many types of charge-transfer processes are possible in open systems. The simplest is between two bound states. 
These are the fundamental processes in many chemical reactions, mixed-valence complexes [38j and biological pro- 
cesses, such as photosynthesis and cell respiration [36|, |37j. A second type of electron transfer occurs between two 
coupled semi- infinite many-electron sub-systems A and B held at different chemical potentials. Direct tunneling takes 
place between two quasi-free states subjected to a thermodynamic driving force stemming from the differences in chem- 
ical potentials. Nonequilibrium single-electron transfer statistics between two leads has been studied extensively over 
the past two decades since it reveals the fundamental quantum effects and direct applications to nanodevices[4l[ |42| . 
Single electron counting has many similarities to the more mature field of single photon counting [43], |44| . The steady 
state current between A and B can be computed using the Landauer-Buttiker (LB) scattering matrix formalism 

[Em 



- dE[f A (E) - f B (E)]J2\Sai(E)\ 



(1) 



where fx(E) — [1 + e /3 (- E_A ' A '-'] _1 is the Fermi function for lead, X = A,B with chemical potential fix, and S a b is a 
scattering matrix element for electron transmission from a -th mode on A to the 6-th mode on lead B. We use h = 1. 

Bardeen's perturbative approach [47| has been very successful for imaging metal surfaces in scanning tunneling 
microscopy (STM). Tursoff and Hamman [48[ had used it for computing the current in STM configurations. Using a 
spherical tip geometry and a generic wavefunction for the metal surface which decays exponentially in the tunneling 
region, they arrived at a simple formula for tunneling current 

^Ka(r ,E F ) (2) 

where a is the density of states (DOS) of the metal surface at the tip position ro and the equilibrium energy Ef- 

Electron transport through a quantum system, e.g. a single molecule or a chain of atoms, connected to two macro- 
scopic leads (A and B) constitutes a third type of charge transfer. In this " molecular- wire" configuration, the electron 
moves sequentially from lead A to the system and then to lead B 1 . In contrast to the tunneling case, where each 
of the two leads is held at equilibrium with its own chemical potential, this case is more complex since it involve a 
non-equilibrium state of the embedded quantum system. 

Electron coupling with molecular vibrations may result in inelastic scattering. These processes contain signatures 
of the bonding environment and the nonequilibrium phonon states of molecules in STM junctions (50j |. Such processes 
are not captured by the simple LB scattering picture. The early NEGF formulation of currents through molecular 
junctions developed by Caroli et al 51] took into account inelastic processes and the nonequilibrium state of the 
quantum system. This approach has since been broadly applied to study the conductance properties of molecules in 
STM [25|,[52,[5E[54[ and molecular wire junctions [55|,[56j], and to the optical properties of semiconductors [j|. Recently 
Cheng et al [391 ] have used the time-dependent density-functional theory to numerically compute the conductance of 
a polyacetylene molecular wire and demonstrated that the NEGF theory works well. 



[1] Note that there is a finite probability for the electron to return to the first lead without tunneling across the junction. Such processes 
are relevant for the current statistics applications at low bias discussed in Section llXI 
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FIG. 1: Time contour for Eq. Physical time (t) runs from bottom to top. Arrows at different times indicate interactions 

of the bra or the ket of the density matrix with an external perturbation, a) Evolution in Hilbert space involves first interacting 
with the ket and then with the bra. Loop time variables (t) are ordered clockwise, n < T2 < T3 < T4 < T5 < re. Interactions are 
time-ordered in physical time within each branch but not between branches (partial time ordering b) In Liouville space, 
both bra and ket evolve forward. All interactions are ordered in physical time (t\ < ti < tz < <4 < ts < te )■ 



Current-carrying molecule may show a negative conductance. The change in resonance condition between two 
conducting orbitals from different parts of a molecule can influence the conductance significantly, which, under certain 
conditions, may induce negative conductance [57], [Elj]. Inelastic effects can also lead to negative conductance. This 
was demonstrated by Galperin et al [5|| using the self-consistent solution for the NEGF in the presence of electron- 
phonon interactions. Recently, Maddox et al [59j have used a many-body expansion of the SNGF to explain the 
hysteresis switching behavior observed [60| in a single magnesium porphine molecule absorbed on NiAl surface in an 
STM junction. In order to analyze inelastic effects in STM imaging of single molecules, Lorente and Persson [53[ 
combined the Tersoff-Hamann theory with the formulation of Caroli et aZ[5l| to compute the fractional change in the 
DOS, Eq. ([2]), perturbatively in the electron-phonon coupling using NEGF. 

The article is organized as follows. The various prescriptions for bookkeeping for time-ordering are first introduced 
without alluding to many-body theory in sections [TT1 and HTT1 The Hilbert-space CTPL technique is presented in Sec. 
HI! and in Sec. IIIII we consider the superoperators and their algebra in Liouville space and introduce the PTBK and 
PTSA prescriptions. In Sec. HVI we treat closed interacting systems using the superoperator approach. This illustrates 
how superoperator time-ordering works. Application is made to Van der Waals' forces. The remainder of the review 
presents applications to open many-body systems. In Sec. [Vj we introduce Fermi superoperators. In Sec. I VII we 
define the many-body SNGF and connect them to retarded, advanced and correlation Green's functions. In Sec. IVII1 
we recast various observables such as the interaction energy, the charge density profile, current induced fluorescence 
and electroluminescence of interacting sub-systems in terms of these Green's functions. In Sec. IVIII1 a closed matrix 
Dyson equation is derived for the SNGF which allows to include electron-electron or electron-phonon interactions 
perturbatively through a self-energy matrix. This is reminiscent of the Martin-Siggia-Rose 61] formulation in classical 
mechanics. Only three Green's functions appear in the PTSA formulation (the fourth vanishes identically). The 
calculation is simplified considerably since the matrix Dyson equation decouples from the outset and the retarded, 
advanced and correlation Green's functions can be calculated separately. In the CTPL approach, in contrast, we 
must first calculate four Keldysh Green's functions and then perform a linear transformation to obtain the PTSA 
in terms of the retarded, advanced and correlation Green's functions. This two-step approach is not required in 
the Liouville-space formulation as the microscopic equations of motion can be constructed directly for the PTSA 
operators. Current-fluctuation statistics in tunneling junctions provides a wealth of information. Extending it to the 
single electron counting regime (as is commonly done for photon counting [43],|44j]) is an important recent development. 
This is connected to the SNGF in Sec IIX1 Superoperator algebra for bosons is introduced in Appendix IA145] . Other 
appendices give technical details. We finally conclude in Sec. |X] 
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II. SCHWINGER'S CLOSED TIME PATH LOOP (CTPL) FORMALISM 

We consider a system described by the Hamiltonian, 

H T (t) =H + V(t), 

where V(t) is an external perturbation that drives the system out of equilibrium, and 

H = H + H' 



(3) 



(4) 



is a sum of a zero order part (Hq) and an interaction part H' which will be treated perturbatively. For t < t , 
V{t) — 0, and the system is described by the canonical density matrix, 



Mo) 



-PH 



We are interested in computing the expectation value of an operator O 

S = Tr[6p(t)], t>t Q . 

We work in the Schrodinger picture where 

p{t) = U{tM)p{h)U ] {tM) 

with the time evolution operator 



U{t,t ) = Texp 



i / drH T (r) 



and its hermitian conjugate 



exp 



i / cLtHt{t) 



(5) 

(6) 
(7) 

(8) 
(9) 



Here T{T^) are the forward (backward) time ordering operators: when acting on a product of operators they rearrange 
them in increasing (decreasing) order of time from right to left. 

We next switch from the Schrodinger to the interaction picture where the time dependence of an operator 0(t) is 
given by 



here 



T {t) = W H (t,t )O{t )U H (t,t ) 



U H (t,t ) = 9(t-t Q )cxp{-iH{t-t Q )} 



(10) 



(11) 



represents the evolution with respect to H rather than Ht- The density- matrix evolution in the interaction picture 
is given by, 



where 



pi(t) = Ui(t,t )p(t )Uj(t,to), 

Ul(t,t ) = fexpj-ijf drVi(r)}. 
Uj{t,t ) = ftexp/z / Vi{t) 



(12) 



to 



are the interaction picture propagators. Equation (| 1 2|) can be recast as 

pi(t) = T c cxp {-if drV/(r)l p(t ) 



c 



(13) 



(14) 
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where Tq denotes time-ordering on the contour C shown in Fig. la. Our observable, Eq. (J6j> , finally becomes 

'fcexpi-i [ drVHT)ld/(t)\. (15) 



I Jc ) I 

Here (• • • ) = Tr[- • • /3(to)] is the average with respect to the density matrix Eq. (J5]) which includes the interactions H' . 
For practical applications, one has to convert this to the expectation value with respect to Hq. The contour in Fig. la 
is then modified to account for the initial correlations [8j, [9j, [62[ which are important only for the transient properties. 
In most physical problems, like non-equilibrium steady state, this modification is not important. Note that the time 
dependence on H appears at three places in Eq. (fT5|) [in Vi(t),Oi(t) and p(to)]- Using a second interaction picture, 
we can switch the time dependence from H to H . The density matrix p(t ) is obtained by adiabatic (Gell-Mann 
law) switching of the interaction starting at t — * — oo, where the system is described by Hq alone. 

p(t ) = Ui(t , - oo)p U J (t , -oo) (16) 

po is the density matrix of the non-interacting system (Eq. |5j) with H replaced by Ho) and 

UiQto,-ao) = Tcxpj-z [° dTH' Ha {r)\ (17) 



H' Ho (r) = Ul o (r,0)H'U Ho (r,0) (18) 
with Ujj B (t,0) = 8(t)e~ lHot is the time evolution with respect to Ho. Equation (|10[) can be similarly expressed as, 

6i(t) = u\{t,to)OH Q {t)tJi(t,to)- (19) 

Substituting Eqs. (fTT>)) and (|19p in (Qjj]) and combining the exponentials (this can be done thanks to the time 
ordering operators which allows us to move operators within a product of operators), we finally obtain 

S = (f c exp |-» J drW Ho (r)| 6 Ho (t)^ . (20) 

where r varies on a contour C which starts at t = — oo and Wf/ (t) = H'jj (t) + Vjj (t). All time dependence 
is now given by Ho and (• • • )o is the trace with respect to po- This formally looks like equilibrium (zero or finite 
temperature) theory. The only difference is that the real time integrals are replaced by contour integrals. Thus in the 
CTPL approach, the non-equilibrium theory is mapped onto the equilibrium one, where standard Feynman diagram 
techniques and Wick's theorem can be applied. 



III. LIOUVILLE-SPACE SUPEROPERATORS AND THEIR ALGEBRA 

The N x N density matrix in Hilbert space is written in Liouville space [H, H3, [fiH as a vector of length TV 2 . With 
any Hilbert space operator A we associate two Liouville space superoperators labeled "left", Al, and "right", Ar. 
These are defined by their actions on any other operator X [23, [28| 

A L X = AX, A R X = XA. (21) 
We also introduce the unitary transformation 

A + = ±={A L +A R ), A_ = ^=(A L -A R ), (22) 

The inverse transformation can be obtained simply by interchanging + and — with L and R, respectively. In the 
following expressions we denote superoperator indices by Greek letters(a, (3, k, rf). These can be either +, — (the PTSA 
representation) or L,R (the PTBK representation). In the +,— representation a single act of A- in Liouville space 
represents the commutation with A in the Hilbert space. Thus the nested commutators that appear in perturbative 
expansions in Hilbert space transform to a compact notation which are more easy to interpret in terms of the double 
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sided Feynman diagrams[20|. Similarly, a single action of A+ in Liouville space corresponds to an anticommutator in 
Hilbcrt space. 

A-X = -j={AX-XA) (23) 
A+X = -L(Ax + xA). (24) 

Note that in the " single time" nomenclature [3, EH one often denotes the left branch by + (positive time for the 
loop) and the right branch by — (negative). Here we denote these by L and R. The "plus" (+) and "minus" (— ) 
are reserved for the sum and difference combinations in the PTSA representation. We believe that this notation does 
most justice to the significance of the two pictures within the density matrix framework. 

With any product of operators in Hilbert space, we can associate corresponding superoperators in Liouville space. 

{AiA r --A k ) L = A iL A jL ■ ■ ■ A kL 

{AiA r --A k ) R = A kR ---A jR A iR . (25) 



This follows directly from the basic definition, Eq. (|2 1 f) . 

The following rules which immediately follow from Eqs. (j22|) and (|25l) allow to convert products of operators directly 
to products of +/— superoperators, 

(AAj)- = — ^ [[A i+ ,A j+ ] + [Ai-,A^] 

+ {A i+ ,A j _} + {A i -,A j+ }} (26) 

{A l A 3 )+ = — = + 

+ [A i+ ,A j .] + [A i .,A j+ ]\ . (27) 

Equations (|2"6"|) and (|27[) may be used to recast functions of Hilbert space operators, such as the Hamiltonian, in terms 
of superoperators in Liouville srjace. 

The time-ordering operator T is a key tool in the Liouville space formalism; when acting on a product of time- 
dependent superoperators, it rearranges them in increasing order of time from right to left. 

fA za (t)A jfj (t') = 6(t' - t)A j0 (t')A ia {t) + 6(t - t')A ia (t)A jP (t') (28) 

Note that, unlike Hilbert space where we need two time-ordering operators to describe the evolution in opposite 
forward (T) and backward (T<) directions (Eqs. (J5J) and ([9])), the Liouville space operator T always acts to its 
right and only forward propagation is required. This facilitates the physical real-time interpretation of all algebraic 
expressions obtained in perturbative expansions. 

Using Eqs. (|2ip and (f2"2")l . it is straightforward to show that 

(A L ) = (A R ) = (A) 

(A + ) = V2(A), (I_)=0. (29) 

where (•) represents a trace which includes the density matrix, (A)=Ti{Ap}. vanishes since it is the tract of a 

commutator. An important consequence of the last equality in Eq. (|2"9")) is that the average of a product of "plus" (+) 
and "minus" (— ) operators with left most "minus" operator vanishes, i.e. (Ax- ■ ■ ■ A na ) = and therefore the average 
of a product of " all-minus" operators vanishes identically (A1-A2- ■ ■ ■ A n J) = 0. The following time-ordered average 
of two superoperators is therefore causal 

(M l+ ft>i^)> = { Jti&m- » t > t m 

For fermion operators A R is defined differently than in Eq. (|2ip in order to maintain simple superoperator commu- 
tation rules. Consequently, some of the above results, in particular Eqs. (|29|) . no longer hold. Fermion superoperators 
will be discussed in Sec. IVl 
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A. The interaction-picture 

We are interested in computing the expectation value, Eq. ([6]), for the system described by the Hamiltonian, Eq. 

The superoperators corresponding to H, Hq,V and VV will be denoted as H,Ho,V and VV, respectively. 
The time-evolution of the density matrix in the Schrodinger picture is given by the Liouville equation 

^p(t) = -i[H T (t),p(t)] (31) 



which in superoperator notation reads, 



The formal solution of Eq. (|32|) is 



where 



dt 



^p(t) = ~iV2H T ^p(t)- (32) 



p(t) = U(t,t )p(t ) (33) 



U(t, t ) = Texp j-n/2 jf drff T _ (r) \ ( ,( /„ ) ( :!4 i 



is the time evolution operator in Liouville space. 

By substituting Eq. © in (J34|), we can switch to the interaction picture 



where 



represents the free evolution, and 



U{t,to) = U (t,to)Ui(t,t ), (35) 
U (t,to) = e(t-t )e- iV§6o -^- to \ (36) 



fexpj-iv^ jT dr>Vi(r)| . (37) 



In the interaction picture, the time dependence of any superoperator is given by, 

Wi(t) = W (t,t )W a (t )tj (t,t ). (38) 

where a = L,R or +, — . Interaction picture superoperators will be denoted by a superscript /, A^. Using this 
notation, our observable, Eq. ([6]), is then given by 

S = Tr[0 I (t)p I (t)} (39) 

where pi(t) — Ui(t,t )p(to) is the density matrix in the interaction picture, pi can be generated from the density 
matrix, p$, of the non-interacting system by switching the interaction W adiabatically starting at to = — oo. 

p I {t) = U I {t,-<x>)p (40) 

Combining Eqs. (j40|) and ([39]) gives 



S = (fcxp |-i\/2 J W^(r)dr| 0j( 



(*) > (41) 



In contrast to Eq. (|20[) , we note that we only need a single time-ordering operator which is always defined in real-time. 
By expanding the exponential, S can be computed perturbatively in the interaction V. 
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IV. SUPEROPERATOR FORMALISM FOR CLOSED INTERACTING SYSTEMS: VAN DER WAALS' 

FORCES 

Here we demonstrate how the superoperator approach is used by applying it to intermolecular forces (29l . l30l . l3ll ] . 
This "single-body" example sets the stage for the many-body applications presented in the remainder of this review. 

Consider two interacting systems of charged particles, A and B, kept at sufficiently large distance such that their 
charges do not overlap and their interaction is purely Coulombic. The total Hamiltonian is 

H = H A + H B + H AB (42) 
where H A H b represent the isolated systems and the coupling is 

H AB = J dv J dv'K{\r-v'\)Q A (v)Q B {v') (43) 

where Qx( r ) is the charge-density operator for system X at point r and KQv — r'|) = l/|r — r'| is the coulomb 
potential. 

The interaction energy is defined as the change in the total energy of the system due to interaction H AB , of the 
system. 



AE = (H ab )ab - (H A ) a - {H B )i 



(44) 



where (• • • )ab is the trace with respect to the ground state density-matrix of the interacting system and (• • • )a is the 
trace with respect to noninteracting density-matrix of A, and similarly for B. 
We introduce a switching parameter < 77 < 1 and define the Hamiltonian 

H V =H A + H B + V H AB . (45) 

Using the Hellman-Feynman theorem, Eq. (|44p gives 







AE = I dT 1 {H AB ) r , = -j=J^ d v (H AB+ ) v (46) 



where (• • • ) v denotes a trace with respect to the density- matrix corresponding to the Hamiltonian H n . The integral 
in Eq. (|46p is due to the adiabatic switching of parameter r\ from zero (noninteracting systems) to one (when 
H ABr) = H AB ). In going from first equality to second we have used Eq. ([2^1) . 

The superoperators, H ABa ,a = +, — corresponding to the coupling, H AB , are obtained using Eqs. (|26[) and Q27[) 
in Eq. (|4*3")l . Since Q A and Q B commute, we get 



H AB — = J dv J dv'K{\v-v'\){Q A {v)Q B {v'))^ 

= -^JdrJ dv'K{\v-v'\)[Q A+ (v)Q B -^) + QA-^)QB + {v') 



(47) 



Similarly, 



H AB+ = J dv J dv'K{\v-Y l \){Q A (r)Q B (v l )) + 

= -L J dv J dv'K(\v-v'\)[Q A+ (v)Q B+ (v')+Q A _(v)Q B _(v')}. (48) 

Substituting Eq. (|48|) in (f44|) and using the fact that (Q A -(v)Q B -{v')) = 0, we obtain 

AE=±j\ v JdvJ dv'KQv - v'\)(Q A+ (v)Q B+ (v')) n . (49) 

We next define the charge-density fluctuation operator corresponding in system X, 5Qx(v,t) = Qx(v,t) — Qx(v), 
where Q(v) is the average charge-density at point r. The total interaction energy in Eq. (|49|) can then be factorized 
into three parts, AE — AE a + AE\ + AE 2 , where 

AE = j dv J dv'K(\v-v'\)Q A (v)Q B (v') (50) 
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is the average (classical) electrostatic interaction energy between A and B while SEi and 6E2 contain the correlated 
density fluctuations. 

d V fdr f dr'Kdr-r'DlQA^iSQB+ir^^ + QB^iSQA+ir))^ 



AEi = 
AE 2 = 



dr, / dr / dr'K(|r-r , |)(5Q A+ (r)(SQ B+ (r'))„ 



(51) 



Using the interaction picture with respect to Hab , each correlation function in Eq. (l51j) can be expressed as 

(SQA+(r)SQ B+ (v')) v 

= (f5& A+ (T,t)8& B+ (T',t)exp{-i^J dT V H AB _(T)yj , (52) 

where the trace is with respect to the direct product density-matrix of systems A and B. By expanding the exponential 
we can factorize the correlation function at each order into a product of density-correlation functions for system A 
and B. 



SQA ai ( r > n) • • ■ SQ Aa „ 0, T n )SQ Ban+1 (r, r n+1 ) ■ ■ ■ 6Q Ban+m (r, t n+n 



(53) 

To lowest order in the interaction K(\r — r'|), AE2 gives the well known McLachlan's expression, AEmc, for the Van 
der Waals' energy pTOj 

AE MC = \ J dr J dr' J dr, J dr\ J dt J dhKQr - r'\)K{\n - r\|) 
x R A ++{r, t; n, t 1 )R B +-(r', t; r\, ti) 

+ Ra+- (r, t; r 1 ,t 1 )R B ++{r' , t; r'i, ti) 

where -R++ and are the correlation and the response functions of density- fluctuations, respectively. 

Rx++{r,t;r',t 1 ) = (f 8Q x +{v, t)SQ x +(r', h)' 



(54) 



x 



Rx+-(^t;r',t 1 ) 



i(f5Qx+(T, t)SQx-(r',h) 



x 



(55) 



Using the Fourier transform f(cu) = J die lut f(t) and the fact that R++ and R^ only depend on the difference of 

their time arguments, Eq. (|54[) can be expressed in the frequency domain as 

AE MC = \ J dr J dr' J dr, J dr\ J ^-K(\r - r'Q^On - r^) 

x R A++ (r,ri;ui)R B+ -(r',r'i;u)) 

+ i?A+-(r,r 1 ;w)i? B++ (r',r' 1 ;w)] . (56) 

This is a general expression (to lowest order in the coupling) for the interaction energy in terms of the charge-de nsity 
fluctuations (R++) and response functions (R^ ) of both systems. Using the fluctuation-dissipation relation[30l l3l| 



R x ++{r,r',u) = 2coth ( ) lmR x +-(r,r' ,uj) 



(57) 



where (3 — l/k B T, the interaction energy can be expressed solely in terms of the density-response R^ (oj) of the two 

systems. 



AEmc = J dr J dr' J dr, J dr\ J ^-K (|r - r'|)A(|n - r^) 

Ra+- (r, r 1 ; w)Imi? B+ _ (r', r\ ; u) 
+ R B +-{r' ,r' 1 ;ui)lmR A +-{r, r^cu) 



x coth 



(3uj 



(58) 
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When the two systems are well separated (compared to their sizes) we can expand the the coupling K |r — r'| in 
multipoles around the position (charge centers) of the two systems and r° B . Hab can be expressed in terms of 
multipoles of the two systems. 

Hab = Jij^ail^bj + K-ijklfJ-aiQbjk ~ QaijUbk] + '•' (59) 

where fi a and Q are dipole and quadruple operators for molecule A and indices k denote the cartesian axis. 

1 



J» = ViV'i- 



r — r' 



1. 



r— r , ,r — r' 

i 



r — r 



(60) 



r— r . .r' — r' 



The leading dipole-dipole term in Eq. (|59|) . which varies as ~ 1/ |r — r'| 3 in Eq. (|59|) . dominates at large separations. 
Assuming for simplicity that the dipoles of systems A and B are aligned along the same cartesian axis, the matrix 
element Jij in Eq. (|60|) reduces to a number equal to 1/R 6 , where R = \r A — r B \ is the distance between the centers 
of the two dipoles. Equation (|56|) now becomes 



1 f duj 

AEjv/c = j^g / — [a A ++(-u)a B +-(u) + a B ++(-u)a A +-(u) 



(61) 



Using the fluctuation-dissipation relation, it can be recast as 



A^mc = -^e J ^ coth (t") [ a A+-(-u)Ima B +-{uj) + a B +- (-u)Ima A+ - (w)] . 



(62) 



and R^ arc now replaced by generalized polarizabilities a++ and a+- 

a x ++(t-t 1 ) = (ffj, x +(t)tix+(ti), 



x 



= a x (t - h) + a x (h - t) 
ax+-(t-h) = -i(fjj,x+(t)(ix-(h) 



x 



6(t-ti)[ax{t-h)-ax(ti-t)} 



(63) 



where a x (t — t±) —: (£i x (t)jl x (ti)) is the dipole correlation function of system X. 

We can expand a ++ and aq in terms of the eigenstates and eigenvalues (\a),u> a ) and (|6), u>b) of systems A and 

B, 



a A ++(t-t r ) = 2j2P(a)\Haa>\ 2 cos(ij aa ,(t-t')) 



a A+ -(t-t r ) = 28(t-t r )J2P(a)\Vaa>\ 2 sm(uaa>(t-t')) 



(64) 
(65) 



where /i aa < = (a\/i\a') and P(a) = e @ Ea . Corresponding expressions for system B can be obtained by replacing 
indices A and a with B and b in Eqs. (f6"4")l and (l6"5")) . This gives, 



a A ++(u) = a' A (uj), lma A +-(u) = «a(w), 
where a' A (u>) and a' A (u>) are the real and imaginary parts of the susceptibility, 

1 1 



a A (w) = P(a) \fi a 



LU - LO aa > +17] LU + LU aa > - if] 



(66) 



(67) 



with q* a (uj) = ola{— <*>). 
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Using Eq. ([67]) , it is straightforward to show that 



a' A (cj) = -(l-e-P») aA {w) (68) 
which gives the fluctuation-dissipation relation 

a' A (u) = coth f^j a A (uj). 

Using Eq. (|69p in Eq. (|64p and taking inverse Fourier transform gives, 



(69) 



a A++ (uj) = 2coth f^T^ Imai + _(w) (70) 

The explicit expressions for q;a++ and a^n are 

a A++ (u>) = '^P{a)\^ aa >\ 2 {8{uj + Lo aa ,) + 5{lo - Lo aa ,)) 

a a' 

a .4 + _H = ^P(a)|/w| 2 \ - 1 —. ■ (71) 

aa 

Since Reoq (to) is symmetric while Imaq (u>) is asymmetric in to, Eq. (|62[) reduces to 

ASmc = ^6 / J coth (^r) Im W+- (w)} ■ (72) 

We can further simplify this expression by noting that Re(a ++ (k>)c*H (w)) is symmetric and contributes to the integral 

only at the pole of coth(/3w/2) for w = and this contribution vanishes if we take the principal part. Thus we get, 



1 f duj ( 8u)\ 
AE MC = j^PP J — coth ( il- 1 aA+-Ma B+ _(u;) 



(73) 



where PP denotes the principal part. Thus interaction energy can be expressed solely in terms of the response 

functions a~\ of the two systems. At high temperatures the integral in Eq. (|73|) can be expanded in terms of 

Matsubara frequencies defied by the poles of coth(/3w/2). This gives 



-. n—oo 



A-Bmc = -p- a A+ -(iu n )a B +-{iuj n ) (74) 



i? 6 

»=0 



where uj n = 2irni/ @ and a ' over the sum indictes a half contribution at the ploe u> = 0. Eq. (I74| is the celebrated 

McLachlan expression [70l | for the interaction energy of two coupled molecules with polarizabilities a A ^ and olb-\ 

In Sec. I VIII we shall follow the same procedure to compute different properties of interacting open fermionic 
many-body systems. 



V. SUPEROPERATOR ALGEBRA FOR FERMION OPERATORS 

To combine the results of sections (jlllj) and (|IV|) with many-body theory, we switch to second quantization. Here we 
consider fermionic systems. Bosonic systems are treated in Appendix |Al Electron creation (annihilation) operators, 
tpi^) satisfy the fermi anti-commutation relations. 

{4>a, 4>t} = <U Wa, = 0, $} = 0. (75) 

The many-body density matrix can be represented in Fock space as p = Y^,mn l-^0(-^l> wnere \M) and \N) are 
many-body basis states with m and n electrons, respectively. These are obtained by multiple operations of the Fermi 
creation operators on the vacuum state 



\M) = ft m ■ ■ ■ $$|0), (M\ = (0\Mj 2 ■ ■ ■ i 



(76) 
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\M) is antisymmetric with respect to the interchange of any two particles, i.e., it must change sign when any two out 
of the m indices are interchanged on the right hand side of Eq. (|76[) . This state can also be expressed as the direct 
product of m single particle states summed over all possible m! permutations 

ml 

|M)=^i 7 J f A|m)...|2)|l) (77) 

l—l 

where |1) = t/>J|Q), etc. are the single particle states and Pj is the permutation operator which generates i successive 
permutations of two particles. 

With each Hilbert space operator ip we associate a pair of superoperators in Liouville space, tpL, and ipR, defined 
through their action on a Liouville space vector [66j 

#l\M)(N\ ee $\M){N\, ft L \M){N\ = ft\M){N\ 
i> R \M){N\ = (-l) m - n \M)(N\i> 



ft R \M){N\ = {-l) m - n+1 \M)(N\ft. (78) 



The definitions for the "right" operators differ from Eq. (]21fl by the (—1) factors. Note that ip R is not the hermitian 
conjugate of ipR as will be shown below. These are introduced in order to account more conveniently for the Fermi 
statistics. These factors are not required for bosons where the many-body state is symmetric with respect to particle 
exchange (Appendix [AJ and we can simply use the single-body definitions, Eq. (|2ip . The superoperators defined in 
Eq. J7ff[ ) satisfy the same anti- commutation relations as their Hilbert space counterparts [Eq. (|75|) ] 

{L,^} = 0, {^ ma ,i>l } = O. (79) 

The fermion superoperators in the PTSA representation are defined by Eqs. (f2"2")l . Substituting (f?9")) in Eqs. 
and (I27| . we obtain for fermion operators (x,( = V>, i/>') 

(xC)+ = ^= [x+C- + x-C+ + (1 - s XtC )] , (xO- = ^= [x+C+ - C-x-] (80) 

where 6 X ^ is unity if \ arL d C are same type (creation or annihilation) operators and zero otherwise. Equations 1791 
hold in both PTBK (a, = L,R) and PTSA (a, /? = +,-) representations. 
From Eqs. ([75]) . we obtain 



= Tr{^}=Tr{^^|M)(iV|} 

MAT 

ee Tr{£(-ir-«|M)TO} 

MiV 

= ^(M + l|M) = Tr{^|M)(iV+l|} 

A/ MN 

= Tr{^ \M)(N\$} = Tr{#} = (81) 



it /AT 



where in going from second to the third line the sum over N can be done using the fact that trace is the sum of the 
diagonal elements of the matrix \M()N + 1|. Using the same arguments we can write 



MN 

ee Tr{^(-l) m -" +1 |Af)(7V|V> t } 

MN 

= Tr{^(-l) m -" +1 |Af)(7V - 1|} 

MN 

= - J2{M-l\M) = -Tr{J2\M){N-l\} 

M MN 

= -Tr{ftp} = -($). (82) 
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Similarly it can be shown that, i^t) — an d {i>\) — 
Using Eqs. flST]) and (gU) we get 



and 



(#-> 



1 

71 
i 

71 

i 

i 

71 



= 
0. 



(83) 



(84) 



Thus ^ satisfies the relations Eqs. (|29p while does not. Note however that any A a of the form of a product of 
Fermi superoperators which contains equal number of V> and xp* always satisfies Eqs. (f!?9")) . since it behaves like a 
boson operator 2 . 

Wick's theorem, Eq. (|B1[) . for products of fermion operators accounts for the antisymmetry of the many-body state 
with respect to the permutation of any two particles (Appendix |B|) . For fermion superoperators this reads, 



(XlaiX2a 2 ' ' ' X(n-l)a„_ 1 Xna„) — ^ (~ 1) P (Xlai X2a 2 ) ' ' 1 (X(n-l)a„_i/" 



(85) 



where the index p denotes the number of permutations required to arrive at the desired pairing. 



VI. SNGF FOR FERMIONS 

Conventional Hilbert-space CTPL is based on four basic closed time path Green's functions [3. FLU [lH, [T5, l68l. [69l| ) . 
The time contour in CTPL is defined by parameterizing the real time in terms of a loop-time parameter r. A two-time 
Green's function, G(ti,T2), is defined on the contour. Depending on which branch of the contour t± and Ti lie, i.e., 
by ordering the two time-points on the contour, one obtains four Green's functions (G > , G < , G T and G T ), each 
describing a distinct physical process. By formulating Green's function theory in Liouville space, we find that G > 
and G < are related to transport while G T and G T represent Fock space coherences [25| . 

We now have all the ingredients necessary to introduce the SNGF 

G^ n (t,t'):=-i{T$ ma (t)ft nf3 (t')), (86) 
where the indices m, n represent the single particle degrees of freedom, and all operators are in Heisenberg picture 

^(^^-^L^- (87) 

The trace (■ ■ •) in Eq. (|86p is with respect to the density matrix of the interacting system (Eq. ([5]) with H given 
by Eq. (|42p). T is a key formal tool in this approach which allows to derive compact expressions for all Green's 
functions and observables. Each permutation of Fermi superoperators, ip a and ^J,, necessary to achieve the desired 
time ordering brings in a factor of (—1). This sign convention associated with time ordering is essential for deriving 
the Dyson equation for the Green's function, Eq. (|56")l . 

In Appendix [Cl we recast these Green's functions in Hilbert space. In Appendix iDl we show that G__, G++ and 
G |- coincide with the retarded, the advanced and the correlation Green's functions, respectively. 



G m I!(M') = -i6{t-t'){{^ m {t)^ n {t')}) (88) 

G™(t,t') = id(t'-t)({^ m (t),^ n (t')}) (89) 

G™(t) = -i{$ m ®JtV)])] (90) 

G![ m = 0. (91) 



[2] This is true for any operator containing product of even number of ordinary fermi operators. 
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Consider two interacting fermion systems A and B described by the Hamiltonian in Eq. ([42]) where the coupling is 
now given by 



Hab = y~] Jabipli'b + h.c. 



ab 



(92) 



The indices a(b) run over the orbital and/or spin degrees of freedom associated with system A(B). Such coupling 
appears in a wide variety of physical systems [671 ]. We now show how properties of the joint system may be expressed 
in terms of one-particle SNGF of the individual systems. 

At t = — oo we assume that the total density matrix is given by a direct product of density matrices of the individual 
systems, p{— oo) — pa® Pb- The density matrix of the interacting system can then be constructed by turning on the 
interaction Hab adiabatically. 



p(t) = U(t, — oo)p(— oo) 



where 



U (t, — oo) = Texp 



-iy/2 



(93) 



(94) 



H AB _(t) is the superoperator corresponding to the Hamiltonian Hab in the interaction representation. 



H 



AB 



_(t) = U\t,0)H AB U(t,0) 



where U o (t,0) — 9(t)e lH o- t and H _ is the superoperator corresponding to Hamiltonian, H 
interaction representation of the density matrix, Eq. (|93p . the SNGF assumes the form, 



C/f (M') = ~i (f^ ma {t)^{t')^ v 



iV2 



drH AB _(r) 



(95) 

H — Hab ■ Using the 
(96) 



where the trace is now taken with respect to non-interacting density matrix p(— oo) = p A ® Pb- 

Equation (|96p is a convenient starting point for computing the Green's functions perturbatively in Hab-- Each 
term can be expressed as a product of SNGF of the individual systems A and B. The Green's functions of individual 
systems can be calculated perturbatively by solving the matrix Dyson equation derived in section [VIIH 



VII. SUPEROPERATOR GREEN'S-FUNCTION EXPRESSIONS FOR OBSERVABLES 

We consider two interacting systems A and B described by the Hamiltonian, Eqs. ([42]) and ([92]) . For observables 
such as the charge density described by operator A which contain an equal number of creation and annihilation 
operators, if) and ijv we have (Al) = (Ar) = (A) and (A-) = 0. Using a perturbative expansion in Hab we can 
express the various observables in terms of Green's functions of the individual systems. 



A. The interaction energy 



The interaction energy AE is defined by Eq. (|4"4")) where Hab is given by Eq. ([9"2"|) . We follow the procedure 
outlined in Sec. HVI to compute AE. We introduce a switching parameter (77) as in Eq. ([45]) . For 77 = 0, H v describes 
the nonintcracting subsystems A and B, whereas for r\ = 1, we get the full Hamiltonian H\ = H . 77 is a convenient 
bookkeeping device. Using the Hclmann-Fcynman theorem, interaction energy can be expressed as in Eq. (|46[) . 

In the interaction representation Eq. J46|) reads 



AE = ^=J^ drj (fH{ B+ (t)exp S-iriV2 J drH AB _(r)^ 



(97) 



By expanding the exponential we can compute the interaction energy perturbatively in the coupling Jab-, Eq. 
To lowest order we get, 



1 



AE = -Im ^ JabJa'. 



'iba'b' 



dio 
2^ 



G°ll(u)G b *' + {u) + G a l%{uj)G w _(uj) 



(98) 
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This extends the McLachlan's expression, Eq. (|6ip . for Van der Waals forces [30ll3lLl7d.l71| to open fcrmionic systems. 

We next assume that system B is infinitely large and therefore unaffected by the coupling with the small system 
A, and remains at thermodynamic equilibrium at all times. This corresponds, e.g. to a molecule (A) chemisorbed on 
a metal surface (B). We further neglect electron-electron (e — e) interactions in system B and set Hb = J2b e b' l Pb^'b- 
Its Green's functions are then given by 



Gf R (t,t') = id w f b e- ie ^-^ 
G b « L {t,t>) = *<M1 - /fOe—^-*') 
G™'_{t,t') = -i5 w 9(t - i') e - <e > <*-*') 

where f b = [l + expjeb — £f}] is the Fermi function and e b is the energy of an orbital (b) of system B. 
Upon substituting in Eq. (|98|) we get 



AE = -Lj2 J a» J *a'b Im I 

aa'b L 



duj G a _l(uj) 

2?7 LJ — £ b + i0~* 



Re[l-2f(e b )]G a + l(e b ) 



(99) 
(100) 
(101) 



(102) 



This expression depends explicitly on the SNGF of the system A. Electron-electron and electron-phonon interactions 
can be included perturbatively using the matrix Dyson equations derived in Sec. I Villi 



B. charge-density redistribution 



We now turn to a different observable, the fluctuations in the charge density Sua of system A due to its interaction 
with B. The local charge density at point r a is given by the diagonal density-matrix element of molecule A, 

n A {r a ) = PA(r a ,r' a )\ K=ra (103) 

where 

p A {v ai v' a ) = {^{r a )^{v' a )) 

= iG LR (v a ,t;v' a ,t). (104) 

is the density-matrix of system A. The charge density fluctuation of system A at point r a is 

5n A (r a ) = 5p{r a ,r' a )\ ra=r , a 

= \PA(r a ,r' a ) - p A(r a ,r' a )\ ra=r , a (105) 

where Sp A is the fluctuation in p A and poa = (^(r a)tp A(r' a )) 'o is the density matrix of the isolated molecule A. 

We shall compute the change in the density matrix of molecule A, 5pA (r a , r ' a ) perturbatively in J a (,. To that end, 
using Eqs. (|96|) and (|104|) . we recast the density matrix of system A in the interaction picture 



PA(r a ,rl) 



U - 



-i {f^ L (r' a ,t)^ R (r a ,t)exp 



-iV2 [ 



drH AB _{r) 



(106) 



By expanding the exponential we can compute pa perturbatively in the coupling between the two systems, Eq. (|92p . 
The first term in the expansion is simply pqa and the second term (which is second order in the coupling) gives the 
fluctuation of the density matrix. 



SpA(r a ,r' a ) 



£ J drJ{ra 



i - r w )J*(r a2 - v b2 ) 



[G_+(w; r' a ,Y al )G ++ (uj; r a2 , r a )G ++ (cj; r&i, r b2 ) 
G-+(uj; r a2 , r a )G — (w; r' a ,r al )G — (w; r 6 i,r 62 ) 
G-+{uj; r bl ,r b2 )G ++ (uj; r a2 , r a )G — (u; r' Q , r a i)] 



(107) 



where dr = dr a idr b idr a 2dr b2 . 
electrons, SNA- 



Integrating Eq. (|107[) over the region A gives the total change in the number of 



SN A = / dr a 5n A (r a ). 



(108) 
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Note that SNa need not be an integer since the two systems can be entangled; an electron can be delocalized across a 
sub-system, giving them a partial charge. Stated differently, the system can have Fock space coherences. Generally, a 
many-body state *ff(N) of the joint system with N electrons can be expressed as combinations of $ a{u) and ^b{N— n) 
as 



{n,p) m (N-n,q) 
n,p,q 



(109) 



where p, q represent various many-body states of A and B with n and N — n electrons. This entangled many-body 
state may not be represented by an ensemble of states with different n, as is commonly done in ensemble DFT 
0, 0,0,0, III- 

As we did for the interaction energy, we can derive a closed expression for the change in the density matrix of 
system A, in the limiting case when system B is large and made of noninteracting electrons. Denoting the single 
particle (orbital) wavefunctions of system B by tp{r b ), we get 



G++(w;r b ,r b ) = 
G_+(w;r b ,r b ) = 

By substituting Eq. pTU|) in ([107)) . we obtain, 



E 



<p b (r b )(p b (r' b ) 
lu — e b — ST] 



i ^2 ( Pb(rb)'fib(r b )[l - 2f b ]8(uj - e b ). 

b 



2^2 I [ dr al dr a2 I b (r al ,r a2 ) 



(l-2/ b )G++( )G—( 

duj ( G- + (oj;r' a ,r al )G++{uj; r Q 2, r a ) 



2tt l lo — e b — ir\ 

G-+(ui; r a2 , r a )G — (u; r' a ,r a2 ) 
(j — eb + if] 



where 



4(r i,r a 2) = J J dr b idr b2 tp b (r b i)tp b (r' b2 )J(r a i - r b i)J*(r a2 - r b2 ). 



(110) 
(111) 



(112) 



(113) 



C. The current in a junction 

The current is defined as the expectation value of the total electron flux from system A to B. The current operator 
is given by as the rate of change of number of electrons in subsystem A. 



(114) 



Substituting Eq. (f9"2")l we obtain 

I = ie \ Jba1pb1pa - Jab1pti>b) ■ (115) 

The current is given by the expectation value of / with respect to the total (interacting) density matrix of the system. 

I(t) = -2eJ ob Im(4^ ) =2eJ ab ReGt R (t,t), (116) 



where in second equality we have used Eq. (|C3[) . 

As explained in Sec. [V] we can expand the Green's function in Eq. (|1 16[) perturbatively in the coupling between 
the two sub-systems. To lowest order this gives 



'iba'b' 



I(t) = e £ J ab r a , b , / dr\Gil(r,t)Gf R (t,T)-G%(r,t)Gi^(t,r) 



(117) 
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When the couplings J a b are real, the current can also be recast in simple form in terms of the PTSA G and G h 
(Eq. (|E30|) V We note that the total current only depends on G LR and G RL \ Gll and G RR do not contribute. This 
can be understood naturally using the density matrix as depicted in the double sided Feynman diagram as shown in 
Fig. 2. Glr : and GflL represent electron transfer between systems (|JV}(JV| changes to |7V+1)(7V+1| or | — \){N— 1|) 
while Gll and G RR only represent Fock space coherences with no change of the number of electrons (|/V)(/V| remains 
unchanged at the end). 



JV-1){JV-1||JV+]}{AM-1 



w 



r 



W 



\n) {n\ \n) (n\ 

G^tj) G^t/) 



|JV-1) 



tut 



n) {n 



t 



t- 

{/V+l| 

w 



\n) {n\ |w) {n\ 

G RR {t,t') 
{f>t) 



W M I A) (N\ 





t 




t 


¥ 


JV+l} 








(N-l 




{ 









n) (n\ \n) (n\ 

Gjtj') Gjt/) 
(t'<t) 



FIG. 2: Double sided Feynman diagrams for the four superoperator Green's functions. \N) represents a many-body state with 
N electrons. Glr and Grl correspond to the transition from \N) to TV — 1) and | TV) to | TV + 1} state, respectively. Gll and 
Grr represent coherence between the two states during time interval t to t' . 



In the long time limit, the system reaches a steady state where the current becomes independent on time and the 
Green's functions only depend on the difference of their time arguments. Using Fourier transformation of G a p(E) = 
J dtG a p{t)e lEt , Eq. (|117|) can be recast in the frequency domain as 



/ / I dE 

aba'b' 



G R a L {E)Gf R {E) - G^ l {E)GVr{E) 



bb' 



(118) 



As done previously, [Eqs. (|M|) - (|10ip . we next assume that sub-system B is a thermodynamically large, free electrons 
system, in this case, the current can be calculated non-perturbatively to all orders in J. This is derived in Appendix 
|E] Equation (|118|) can then be replaced by Eq. (|E25p 



dE 



m a R {E)G R a L (E) - V R \{E)Gf R {E) 



(119) 



where £J"g ,a, f3 = L, R are the elements of the self-energy matrix for sub-system A due to interaction with B. 



?>Tr{E) = ir±,f B (E), = -ir£,[l - f B (E)} 

^Tl(E) = - l -Tt,[l - Vb{E% Tfj^iE) = |r^,[l - 2f B {E)] 



(120) 
(121) 



with r^ a , = 2tt ^2 bb , J a bJa'b' a B: where gb is the density of states (assumed to be independent of energy) of sub-system 
B. Unlike Eq. (JTT5J) , the Green's functions in Eq. (|119|) now contain coupling to the leads to all orders through the 
self-energies. These can be evaluated by solving the matrix Dyson equation as given in Sec. IVIIII Under certain 
approximations, Eq. (fTT9")) reduces to the LB form, Eq. Q. This is shown in Eqs. (|E32)) - (|E33)l . 



D. Fluorescence induced by an electron or a hole transfer 

Consider an electron transferred to a molecule attached to two leads (a molecular wire) or by an STM probe by 
applying an external potential. This electron may decay radiatively to one of the lower energy states by emitting 
a photon, before finally exiting the molecule. We denote this optical signal current-induced-fluorescence (CIF). The 
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reverse process can occur when an electron leaves the molecule, creating a hole. An electron in one of the higher lying 
bound states can decay radiatively to combine with this hole. Thus the CIF signal can come either from a positively 
or from a negatively charged molecule where the charge is created by the interaction with external leads. This is 
different from ordinary (laser-induced) fluorescence (LIF) [201 where an electron-hole pair, created by the interaction 
with the laser field, recombincs to produce a photon and there is no electron transfer from/to the molecule. Thus 
CIF involves many-body states of the molecule with different total charges whereas LIF only depends on states with 
the same charge. 

The system is described by the Hamiltonian H = H m + Hf + Ha + Hb + H m j + H' where H m , Hf, Ha and Hb 
represent the independent molecule, radiation field, lead A and lead B, respectively. H m f is the coupling between the 
radiation field and the molecule 

H m f =a\B + B^a s . (122) 

a\(a s ) are, respectively, the creation (annihilation) operators for the s-th mode of the scattered field and W(B) is an 
exciton (electron-hole pair) operator for the molecule 

& = J2^Hi, B = J2»jii>li>i- (123) 

fXij is the dipole matrix element between the orbital i and the orbital j. 
The molecule-lead coupling H 1 is 

h'= J2 J ^i^ +h - c - ( 124 ) 

xeA.B i 

The signal is defined by the expectation value of the flux operator for the scattered photon modes, 

S(u s ,t) = (j t ata s ) =i([H,aia s }). (125) 

The angular brackets denote the trace with respect to the full density matrix of the field+molecule+leads system 

Expanding Eq. (|125p perturbatively in the emitted field and in the lead-molecule coupling, the CIF signal from a 
negatively charged molecule can be expressed in terms of Liouville space Green's functions. Depending on the time 
ordering of various superoperators corresponding to B, a and if), the density matrix evolves through different LSPs. 
For the signal to be finite, we must first prepare the molecule in one of its excited states by transferring either an 




FIG. 3: Three Liouville space pathways that contribute to the fluorescence from a negatively charged molecule, Eq. (1126[1 . 
States \a), \b), |c) and \d) denote different charge states of the molecule. Time increases from bottom to top. 

electron or a hole, before bringing it down to lower excited states. This results in the cancelation of many terms in 
the perturbative expansion of Eq. (|125[) . Only three LSPs (plus their complex conjugates) contribute to the signal. 
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These can be expressed by the double sided Feynman diagrams shown in Fig. 3. The fluorescence from a negatively 
charged molecule is finally given by [26[ 



^(w.) = *e]T J2 J2a(E F + eV)J k {E F + eV)J l (E F + eV) 

i<j i' <j' kl 



X jl^j j-Li' j' 



(It 



T,t 



d T . e -iu>,T e -i(EF+eV){n.-Tx) 



if 



i' k 



r) | G% L (t - r, r x )Gl L {r 2 , t) + G» R (t - r, t)G^(r 2 ,n) 



- G^(t,t-r)( 
+ G RL (T 2 ,t - t) (&/ R (t - r, t)G R k L (t, n) + - r, 7i)G^(t, t+)) 



(126) 



where cr(-E) = — ImG--(E)/iv is the density of states of the lead from which electron has been transferred to the 
molecule, evaluated at energy Ep + eV. Jk{E) is the coupling of the kth molecular orbital with the state of the lead 
at energy E. A similar expression can be obtained for the signal from a positively charged molecule (26l|. 



E. Electroluminescence: Fluorescence induced by electron and hole transfer 

Electroluminescence is the recombination of an electron and a hole, injected separately by external charge sources. 
This is a higher order process than described by Eq. (|126p . as both electron and hole are injected into the molecule 
from the leads. In a molecule-lead configuration, as discussed in Eqs. (|122[) and (| 1 24[> . a lowest order perturbative 
expression can be obtained by expanding Eq. () 1 25(1 to fourth order in the left and/or right lead couplings. This 
contains a large number of terms. The process is simpler in the high bias limit (V > fc^T) where an electron is 
transferred from the left lead to the molecule while a hole is transferred from the right lead simultaneously and the 
reverse processes are not allowed. We then get 

S EL (co s ,t) = -8Rc J dr J dn J dT 2 J dT 3 J dT 4 J2 /«(! - h)e- l£a(Tl - T2) 

x e ie b {T S -T A ) e -i^{t-r) J ia J lb J aj J bk 
ijkl 

x (f^M^W^Jn)^^)^^)^) (127) 

where the trace is with respect to the density matrix of the free molecule. The correlation function inside the integral 
can be computed using the many-body states of the molecule. For a free-electron model it can be expressed as sum 
of the products of the Green's functions for the molecule using Wick's theorem. 



VIII. DYSON EQUATIONS FOR THE RETARDED, ADVANCED AND CORRELATION GREEN'S 

FUNCTIONS 



In the previous section we demonstrated how various observables of the coupled system can be expressed in terms 
of the Green's functions of the individual sub-systems, Gap. Here we show how to compute these Green's functions 
with many-body (e.g. electron-electron and electron-phonon) interactions. Using the superoperator approach, a 
closed equation of motion (EOM) for the Green's functions is derived directly in PTSA representation. The retarded 
Green's function satisfies an independent equation (decoupled from the other Green's functions). Nonlinear effects 
due to interactions with the environment arc included in the self-e nerg y which can be computed perturbatively in the 
interaction picture or using the generating function technique [2lL l83l | 

The Hamiltonian of each closed system (A or B) will be partitioned as 

H = H a + H int (128) 

where 



X 



(129) 
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is the free-electron part where x denotes the electronic degrees of freedom, e.g., orbital, spin, wavevector,etc. Hint 
represents the many-body interactions in the system. We need not specify the form of H int at this point. 
In terms of the superoperator PTSA representation, the total Liouville operator is 



where L corresponds to H 



and 



L = L Q + Li nt (130) 



l o = ^2 ex ($1+$*+ ~ ( l31 ) 



V2H int -. (132) 



We shall construct the equation of motion (EOM) for the superoperator, %pi a starting with the superoperator 
Heisenberg equation 

d m 

-Q^ia = i[L,ipxa] (133) 

where [L,ip X a] = Ltp xa — ip xa L. Substituting Eq. (|130| in (|133p and using the commutation rules Eq. ()79|1 of Fermi 
superoperators, we get 

= -i£xi>xa - i\2[H- , ^ xa \ . (134) 
Differentiating Eq. (fB"6")) and using (|134p . we obtain the EOM for the Green's functions, 

(i^ t -e x \Gli{t,t')=5 xx ,6 a p5{t-t')-iV2{T^^ (135) 

When the many-body interaction Hint is turned off, Eq. (|135| reduces to 

(i^ - e x } G° a f\t,t') = <W<W<5(*-*')- (136) 

Equation (|136[) defines the noninteracting Green's functions, G° a g. 

Equation ()135() can be recast in the Dyson form (repeated indices are summed over). 

G(t, = G°(i, t') + G°(t, t'OSCt'r, t' 2 )G(t' 2 , t') (137) 

where G is the 2x2 matrix of Greens functions, G a/3 and the self-energy matrix S is defined by the equation 

Ki>(^'i)G% x \t\x) = -tV2(f[H^t)j XQ mi (t')) 

= -^/f[^(t),^Jt)]^(*Oexp|-^£dr#£(r)|\ 

(138) 

(■ ■ • )o is the trace with respect to the non- interacting density matrix of the sub-system and the time dependence is 
in the interaction picture with respect to H . 

The Dyson equation [Eq. (|137p ] can be written explicitly in the matrix form with 



where we have used the fact that G + = 0, Eq. (1911) . 
The bottom left matrix element in Eq. (|139|) gives 



G° £ + _G__ = 0. (140) 
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Since G++ and G__ do not generally vanish, Eq. (] 140[) implies that S_| = 0. Equation (I139|) then gives the retarded 

G Green's functions thus satisfy their own Dyson equations. 

G__ = G°__ + G°__£__ G__ . (141) 

Similarly we get for the advanced Green's function 

G ++ = G a ++ + G° ++ Z ++ G ++ . (142) 

The correlation function G can then be expressed in terms of the retarded and the advanced Green's functions and 
their self-energies 

G_+ = G°_ + + G°_ + £ ++ G ++ + G°_ _ (S— G_+ + E-+G++) . (143) 

Keldysh had derived Eqs. (|141|) - (| 143[) [2|,|8| by starting with the coupled equations for the four Hilbert space Green's 
functions in the PTBK representation (which in superoperator notation correspond to Gll, Grr, Grl and Glr) and 
then decouple them by transforming to the PTSA representation. The superoperator approach allows us to work with 
the PTSA representation from the outset. The derivation is much more compact and physically-transparent. 

For completeness, we present the corresponding equations for bosons in Appendix [A"l 



IX. CURRENT-FLUCTUATION STATISTICS 



So far we have discussed average observables, such as the total energy, change in the charge density and the current 
in interacting open systems. Fluctuations may be measured in electron tunneling between two conductors and electron 
transport through quantum junctions Much theoretical [zH and experimental [79j effort has been devoted recently 
to the statistical behavior of the transport through such systems. Here we do not review these works but rather give 
a flavor of it. More detailed derivations are given in Ref. [8Cl l8ll|. 



A. Electron-tunneling statistics 

We again consider the two interacting systems A and B described by the Hamiltonian, Eqs. (|42| and (|92|) . The 
number k of electrons transferred in time t can be computed by measuring the number of electrons in one of the 
systems at time t = and again at time t. The difference of the two measurements, k, is a fluctuating variable which 
changes as the measurement is repeated over the same time interval. Thus k has a distribution P{k, t) whose first 
moment gives the average current, Eq. (I117j) . For small applied external bias, the charge transfer can reverse its 
course and these fluctuations may become non-negligible. 

We shall compute the statistics by using the generating function associated with P(k, t) 

x (A,i)^5>* Afc P(M). (144) 

k 

The GF is a continuous function of the parameter A, which is conjugate to the net number of electron transfers 
k. It is convenient for computing moments of various orders by taking derivatives of GF with respect to A. Working 
directly with P(k,t) requires summations over the discrete variable k from — oo to oo. P(k,t) is obtained by the 
inverse Fourier transform of the generating function 

/•27T J\ 

P(k,t) = J -e-^ k X {\t). (145) 

The n-th moment of P(k,t) is given by 

(146) 

The n'th order cumulants are certain combinations of the moments of order n and lower. The first cumulant C^ 1 ' = 
(k(t)) is the first moment and gives the average number of particles transferred from A to B. The second cumulant 
C( 2 ) = ((k(t) - (k(t))) 2 ) is the simplest measure of fluctuations of k around its average value. The third cumulant 
C( 3 ) = ((k(t) — (k(t))) 3 ) gives the skewness (asymmetry) of P(k,t). For a Gaussian distribution all cumulates higher 



(*"(*)> = H) n «^x(A,f) 
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than the first two vanish. The higher cumulants thus measure the non Gaussian nature of the distribution. All 
cumulants of the current statistics can be computed directly by taking the derivatives of the cumulant generating 
function, Z(X,t) = log%(A,i), with respect to A. The n-th cumulant is given by 

(147) 

A=0 

For our model, the generating function, Eq. (1144[) . can be expressed as[8l| 

x(A,t)= r^V(^) (148) 
Jo 2lT 

with 

Z(X, A; t) = log [Tr {Qj(X, A, t)p }] (149) 

and Qi is the interaction picture propagator given in Appendix [Fl 

We shall evaluate Eq. (|149|) perturbatively to second order in the coupling, Eq. (f92|) . As was done in Sec. IVII in 
the infinite past, the density matrix of the two systems is factorized, po — Pa® Pb- The interactions are then built-in 
by adiabatic switching. We then express the cumulant GF, Eq. (|149[) . in terms of the Green's functions of systems A 
and B (see Appendix |F|). 

Z{\,A;t) = 2zsinQ^ J dr x dr 2 {0(n).F(A, n, r 2 ) 

+ S(-T 1 )[e- iA f(0,T ll T 2 )+e iA f(0,T 2 ,r 1 )]} (150) 

with 

T{X,t,t') = ^2 JabJa'b' 
aa'bb' 

All Green's functions appearing in Eqs. (|107|) and (|151[) correspond to the isolated sub-systems. These can 

include many-body interactions e — e or e — p. The Green's functions then include self-energy due to many-body 
interactions as shown in the previous section. The self-energy to first order in e — e interactions is given in Appendix 

K3 



CW(t) = (-if-2(A,t) 



5 ^(M')^(t',t) 



(151) 



B. Current fluctuations: electron-transport statistics 

We now consider electron transport between two leads A and B through a quantum system such as a molecule or 
a quantum dot. The total system Hamiltonian is 

H t = Ha + Hb+Hs + V (152) 

where H a and Hb represent Hamiltonian of two leads A and B, and H$ is the system Hamiltonian. We assume a 
free- electron gas model for all sub-systems. 

H x = e * X = A > B > S - ( 153 ) 

xex 

V represents the coupling Hamiltonian between the molecule and the leads. 

V= J2 JxM4>i + h.c. (154) 

xeA.B 

where J X i is the coupling between the lead (x) and the system (i) orbitals. 

A simple perturbative calculation in the coupling, as was done in going from Eq. (1149)) to (|F5[) . will miss the 
non-equilibrium state of the embedded system. However, using the path- integral technique|81(, the cumulant GF 
Z{X, A) can be expressed in terms of SNGF of the molecule. The final expression for the GF involves determinant of 
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the matrix Green's functions renormalized due to the coupling with A and B. For long counting times, the cumulant 
GF is given by [U 

Z(X, A) = J du lnDet[G _1 (w,A)] (155) 

where the 2x2 matrix G(u>, A) is obtained by solving a Dyson matrix like equation [Eq. (|139[) ] with the following 
self-energy matrix 

+ '-Ti, [e iX (l - f A {u>)) + e-^f A {Lo)} (156) 

-^u'-^ti' [e iA (l-/„M)+e-^/AM] (157) 

"5 r «' " " - ( e_<A - !)/aH)] ( 158 ) 

- i r«(2/ fl ( £( ;)-i)+ir < 4, 

[(e lA + 1)(1 - - (e~ iA + 1)/aM] (159) 

here T^, = 27r^ I , GX J X iJi> x 5(oj — e x ) comes from interactions with the leads. When A = 0, = as discussed 

previously (Eq. (|140[) ). and £__, £++ and S |_ reduce to ordinary retarded, advanced and correlation self-energies 

due to leads-molecule interaction, respectively. 

X. CONCLUSIONS 

We have presented a superoperator Liouville space nonequilibrium Green's function (SNGF) theory in the PTSA 
(+/— ) representation for interacting open many-electron systems. 

By constructing the microscopic equations of motion for superoperators in PTSA representation, the response 
and correlation functions are treated along the same footing in Liouville space. This is reminiscent of the MSR 
formulation [6l| of classical statistical mechanics. The various Green's functions of the CTPL formulation appear 
naturally as a consequence of a simple superoperator time ordering operation which controls the evolution of the 
bra and the ket of the density matrix in Liouville space. The interaction energy and the change in density matrix 
of individual sub-systems are directly recast in terms of the retarded, advanced and correlation Green's functions. 
Electron-transfer statistics in both tunneling and transport may be formulated in terms of the SNGF. Many-body 
effects (such as e — e interactions) can be included perturbatively through self-energies by solving the Dyson equations. 

A notable advantage of the PTSA representation is that since = 0, the matrix Dyson equations of the other three 

Green's functions are decoupled. Dyson equations for the retarded the and advanced Green's functions are derived 
without computing first the matrix Dyson equation for four Keldysh Green's functions in PTBK representation. This 
is because the linear transformation which results in PTSA operators can be performed at the operator level rather 
than on the Green's functions or the generating function. In Liouville space we can work directly with the operators in 
PTSA representation. In Hilbert space (or PTBK), in contrast, one has to first compute the four (lesser, greater, time 
ordered, anti-time ordered) Green's functions and only then obtain the retarded and advanced Greens functions using 
the matrix transformation. This is not necessary in Liouville space. The superoperator Green's function formulation 
for bosons is given in Appendix [A"l for completeness. 
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APPENDIX A: SNGF FOR BOSONS 

We denote the boson creation and annihilation operators by 4>i and <j>\ , respectively. These satisfy the commutation 
relations 

= Sij, [hh\ = [#,#] = (Al) 



£«' + (u,,A) = 
£«'_(u;,A) = 
S+'-KA) = 
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the corresponding superoperators, <fr a and c/y a , in Liouville space are denned by Eq. ([21 



x = <j>x, 4>[x = 4>x 



(A2) 



A notable difference between Eqs. (|A2j> and (|78[) is the absence of the (—1) factors for bosons. Boson superoperators 
in PTSA representation ("+" and "-") are defined by Eqs. ([22]) and satisfy commutation relations similar to Eqs. 



[4>ma,4> n fi] = S mn S a p, [4> ma , (f> n p] = [^ma'^nffl = °' 

Using Eqs. (|A3j in (26]) and (H}, we get 

1 



j S>- = T2 
f+, 1 



n/2 



(A3) 

(A4) 
(A5) 



The corresponding expressions for (<f>i<j)j)_ and {4>i4>j)+ can be obtained from Eq. (j A4[l and (|A5|) by simply neglecting 
the <5y terms inside the brackets. 

In analogy to Eq. (|86|). the boson superoperator Green's functions are defined as, 

D^(t,t')=-i(f} ma (t)ft nl3 (t')). (A6) 

Using Eqs. (|A2|) and (|A6|) . it is easy to show that , D ++ and D j_ represent the retarded, advanced and 



correlation Green's functions, and D- 



0. This gives the boson analogue of Eqs. 



Df_(t, t') 
D l l+{t, t') 
D l l + (t,t>) = - 
D l l_(t,t') = 



-*0(i-O<[&(t),$(f)]) 
i6{t'-t){[Ut),4>){t')]) 

-i{{Ut)Mf)}) 



(A7) 
(A8) 
(A9) 
(A10) 



Unlike Fcrmions, in this case = while D_| represents the retarded Green's function. This comes from the 

different commutation relations, Eq. (|A3|) and (|75|) . The different particle-statistics for bosons affects the form of the 
Dyson equation (Eq. (|139[) ). 

The Dyson equation for bosons, D = D° + D°H D, can be derived following the same steps which lead to Eq. 
(p9]l . Here 



D 





D+- 



D° 







D° ++ 



'--+ 



(All) 



where D° is the Green's function for the reference system and S is the self-energy due to many-body interactions. 

Comparing the top right matrix elements on both sides then gives, S ++ = 0. We finally get the boson analogues of 
Eqs. ipo] ) -(|T33" l) . 



D+- = Dl 
D, = D°__ 



D++ = D\ + ~ + -D_ + + D, 



(A12) 
(A13) 
(A14) 



These equations are decoupled and each of the Green's functions can be computed separately. and D |_ represent 

the retarded and the advanced boson Green's functions. 



APPENDIX B: WICK'S THEOREM 



Wick's theorem states that the expectation value of the product of n fermion field operators with respect to a 
many-body state, represented by a single Slater determinant, can be expressed as the sum of the expectation values 
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of products of all possible pairs of operators @, [H, H2, 83] . This theorem provides a practical tool for factorizing the 
expectation value of a time ordered product of operators to simpler quantities. Note that expectation value vanishes 
unless n is even. 

This theorem can be derived as follows. Since a single determinant represents the eigenstate of corresponds to a 
noninteracting Hamiltonian, Hq = ■ KijAiAj, the expectation of a product of n operators (AiA 2 ■ ■ ■ A n ) involves 

an integral weighted by a Gaussian operator ~ exp ^— (3 AiKijAA . Using properties of Gaussian integrals, it can 
be factorized as [l|| HI] . 

(A X A 2 ■ ■ ■ A n ) = ^(A 1 A 2 ) ■ ■ ■ (i„_ri n ) (Bl) 
p 

where the sum over P runs over all possible pairings of {1, 2, • • • , n}, Equation (|B1[) remains valid also for superop- 
erators. By simply replacing all A n >s by A nai a — L,R, +, — . 

(A la A 2 p ■ ■ ■ i„ 7 ) = ^2(A la A 20 ) ■ ■ ■ (i„_i 5 I„ 7 ). (B2) 
p 

Since each permutation of fermion superoperators brings a minus sign, the Wick's theorem for fermion superoperators 
takes the form, 

(A la A 20 ■ ■ ■ A m ) = ^2(-iy{Ai a A 2p ) ■ ■ ■ (A n _ ls A ni ) (B3) 
p 

where p is the number of permutations required to get the right pairing. The same theorem applies to bosons by 
simply eliminating the (— l) p factors. 



APPENDIX C: CONNECTING THE CLOSED-TIME LOOP WITH THE REAL TIME GREEN'S 

FUNCTIONS 



Standard CTPL Green's function theory is formulated in terms of four Hilbert space Green functions: time ordered 
(G T ), anti-time ordered (G T ), greater (G > ) and lesser (G < )[3, 0, 0. These are defined in the Heisenberg picture as, 

-i(TV)(t 1 )^t( t2 )) 

-imi)^(t 2 )) 



G T (h,t 2 ) 

G f (h,t 2 ) 

G > (t 1 ,t 2 ) 
G<(h,t 2 ) 



9(t 2 
0(h 



(CI) 



T (T) is the Hilbert space time (anti-time) ordering operator, Eqs. ([8]) and ([9]). 

The four Liouville space Green's functions which appear naturally in the superoperator formulation are 

G LL (t u t 2 ) = -»C%(ti)$,(t 2 )}, G RR {t l ,t 2 ) = -i{fj> R (t l )i> ] R (t 2 )) 



G LR (h,t 2 ) = -i{f^ L {tx)i) R {t 2 )), G RL {t u t 2 ) = -i{Ti>R{tiW L {t 2 )) (C2) 
To establish the connection between Eqs. (|C1[) and (|C2[) . we shall convert the superoperators back to ordinary 



operators by using their definitions and anti-commutation relations. For Glr and G R l, we obtain, 

G LR (t 1 ,t 2 ) = -iTr{fif, L (ti)$ R {t 2 )p} 

= -i [0(h - t 2 )Tr{Mti)^ R (t2)p} - 0(t 2 - ii)T*{$k(t2)$L(tiM 

= iTr{^(t x )(4^t 2 )} = i{ft{tS{t x )) 

= G<(h,t 2 ) 

G RL (h,t 2 ) = -iTrif^RitSUt^p} 

= -i [e{h - t 2 )Tr{i> R (h)$l(t 2 )p} - 6{t 2 

= iTr{ft{t 2 )p4>(h)} = imi)^{t 2 )) 

= -G>(h,t 2 ) 



(C3) 



t^^Ut^R^p} 



(C4) 
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where p is the fully interacting many-body density matrix. 

For Gll and G RR we need to distinguish between two cases, 
(i). For fx > t 2 , we get, 



G LL {h,t 2 ) = -iTr{T$ L (ii)$,(*2M 

= -iTr{i>(ti)ft(t 2 )p} = -t$(fi)$t(t 2 )) 

= -.^(^(t!)} = -i(v;t( t2 )^( tl )) 



(ii) For the opposite case, ti < i 2 , we get, 



(C5) 



G LL {ti,t 2 ) = -iTr{f$ L (ti)ft L (t 2 )p} 

= m{^(t 2 )ip(ti)p} = (t 2 )^(h)) 

G RR (h,t 2 ) = -zTr{fy^(fi)V^(i 2 M 

= iTr{(4(ti)ft{tz)} = *<#i)$ + (*2)) 



Combining Eqs. (|C5|) and (|C6[) we can write, 

G LL {tiM) = ^Tr{T^ L (fi)^(t 2 )p} 

= G T (h,t 2 ) 
G RR {t x ,t 2 ) = -iTr{fi, R (tS R (t2)p} 

= -i [e(h - t 2 )(^(t 2 )^(h)) - 9(t 2 - ti)$(ti)ft(t 2 )) 

= -G f (t u t 2 ) 

Equations (|C3[) - (|C7[) show the equivalence of the Hilbert and the Liouville space Green's functions. 



(C6) 



(C7) 



APPENDIX D: THE RETARDED, ADVANCED AND CORRELATION LIOUVILLE-SPACE GREEN'S 

FUNCTIONS 



Here we show that the SNGFs G , G ++ and G h defined as, 

G%(t,t') = -i{T$Ut)i>Ut')) 

coincide with the retarded, advanced and correlation functions, respectively. 
Starting with G , we expand the time ordering operator as 

G l l_(t,t') = -t(f^(t)ty_(t>)) 

0(t - t')$i- (t)$_ (0) - 0(t> - t) (t'^i- (t)) 



Substituting -0- = (V>l — ip R )/V2, we can write 



(f)) = 5 [(fc(i)-^(t))ft(tO-ft(0))] 



{^L{t)^) R {t'))-{^ R {t)^] L {t')). 



Using the definitions of -ipL and ijj R [Eq. ([78)) ] . we get 



(Dl) 



(D2) 



(D3) 



(D4) 



Similarly, 



= (i(i')fe(')) + $U(t)$iR(t)) 

- (^(O^L(«))-(^(O^Ii(«)) 
EE <^(^(i)) + (^j(0) 



Substituting from Eqs. (|D4|) and (|D5|l in (|D2|) . we get Eq 
Similarly, 



for the retarded Green's function. 



Since 



and 



6(t - t'){^ + (t)^ + (t'))9{t' - t)(i>Ut')i> i+ (t)) 



= o 



<{&(*), #(f)}>, 



substituting these in Eq. (|D6|) . we get Eq. 
We next turn to G u 



for the advanced Greens function. 



G l -+{t, t') = -i<f&_(t)$ > 



0(i - o<&- (*)$+(*')) - - 



Substituting for tp_ , ip + in terms of ^ and ip R gives 

= \{{^L{t)-^ R m^ L (t')+^R{t'))) 

X^dt)^')) + {^L{t)^] R {t')) 



(i> iR (t)ft jL (t')) - (&B(t)$ fl (f)> 



(W^(0)-W}(f*(t)) 



t/v\ 
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Similarly we can simplify the other correlation function 



■;<.'■;+(/')<''; (/.)) = + v-h(O)^/-(0 - tHu(t))) 

Ir.vt^L- ,.y.t 



t ^ 



~{[UtUW)\). 



(Dll) 



Equation (|90)) is obtained by substituting Eqs. (|E)1 1|) and (|D10[) into (|D9|) . This function is non-causal. 
Finally we consider 

= -% [e(t - t')$ i+ {t)i>]_{t')) - 6(t' - t)(i>t_(t')i> i+ (t)) 

Since 



(D12) 



2 

l r 



+ fe(*)C(*')> - (vW^L(O) 



(iM*)$(0> + W(0^(t)> 



1 

2 

(^(t)$(i0)-<$(o^)) 



(D13) 



Similarly it can be shown that (V>j_ (*)$£+(£')) = 0. Equation (f9Tj) follows from Eq. (|D12|) . 

APPENDIX E: CURRENTS THROUGH MOLECULAR JUNCTIONS WITH FREE-ELECTRON LEADS 



In Sec. I VII CI we assumed a small quantum system A coupled to an infinite system B. This corresponds to molecular 
wire or STM experiments where a single molecule (or a quantum dot), the "system", is coupled to two metal leads 
which are large compared to the system. In this appendix we show how by adopting a free electron model of the 
leads, the current can be calculated nonperturbatively in the coupling with the leads. The total Hamiltonian is 



Ht = H$ + H / 



V. 



(El) 



Ha and Hb are the free leads Hamiltonians given in Eq. (|153p . V is the coupling of the system with the leads, Eq. 
(|154[) . The system Hamiltonian in this case is quite general and can include electron-electron and/or electron- phonon 
interactions. 



ijkl 



(E2) 



where the fermion operators i/>^(V>) satisfy the anti-commutation relations, Eq. 1751 Indices (i,j,k,l) and (m,n) 
represent the system electronic states and system phonon states respectively, e, and fl m are corresponding energies. 
Xij and Vijki are the e — p and e — e interaction matrix elements, ^((^m) are phonon creation (annihilation) operator 
for the m-th normal mode, with frequency O m . These operators 4> satisfy the boson commutation relations, Eqs. 
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The current from lead A (left) to the system is defined as the rate of change of the total charge of the lead, 

J W = E |<^> =i^2(W T ,N A (t)]) (E3) 
i i 

where N A = e J2 a V'lV'a is the number operator for the left lead. We substitute the Hamiltonian and calculate the 
commutator. Na commutes with all terms in the Hamiltonian, except the interaction V, Eq. (] 154|) . We obtain 

I(t) = -2elm 1^ J W (V4 (t)^«(t))r| (E4) 

where (• • • )t is the trace with respect to the total (system+leads) density matrix. A Similar expression can be written 
for the current from the right lead B to the system by changing a with b in Eq. (|E4p . At steady state the two currents 
are same. 

The Liouville space Green's functions are connected with their Hilbert space counterparts in Appendix [U] The 
current, Eq. (|E4|) , can be expressed in terms of the Liouville space Greens function as 

I(t) = 2eRe{J ia (t)Gf R (t, t)}. (E5) 

In Eq. (|E5[) . the Green's function GT™ is defined in the joint leads+system space. We compute GLFt(t,t') and then 
set t = t' . In a recent work Rabani [85| has used a numerical path- integral approach to directly compute the Green's 
function Glr (hence the current). 

We shall compute this Green's function using the equation of motion technique. To that end, we construct the 
equation of motion for ipaa(t)(a — L, R) using the Heisenberg equation for superoperators 

^j$aa(t) =lV2[H_J aa (t)] (E6) 

where on the right-hand side 

[fl_,$ oa (t)] = fl-$aa(t) - i>aa{t)H- (E7) 

and H- = (H L - H R )/y/2 can be expressed in terms of ipL and ipR using Eq. (|25|) . Equation (|E6j) then gives 

»^a«(t) = ^ aa {t) + J* at {t)^ a (t). (E8) 

i 

Taking the time derivative of Eq. (|86[) with x = l,y = i and fi = a,v = /3, we get 

8 



%{t,t') = S(t - t')5 a p5 ai + (T (J- t $a a (t)\ $ (t')). (E9) 



The 5(t — t') factor comes from the time derivation of the step function (which originates from T) and the Kronecker 
delta functions result from the commutation relations, Eqs. (|T9"|) . Using Eq. (|E8|) in (|E9p we get 

(i§- t - ta)G%(t, = E W) G Zf>& ( E1 °) 
j 

Note that a ^ i, since they belong to different regions of the total system. Thus the first term on the r.h.s. in Eq. 
(|E9p docs not contribute in this case. To zeroth order in lead-system interaction, the Green's function for the leads 
is defined as, 

{ i ~ £a)5 ^' ( *' <,) = ^-*')<w<w. (En) 

Substituting this in Eq. ijElO]) . we can write 

G&(*,0 = gT^h^A^Gf^t'). (E12) 
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For Glr this gives, 



The current, Eq. (|E5[) . then becomes 



I(t) = 2eRe £ J d^J^g^i^G^^t). 

ii' aa' (3^ 



We define the Fourier transform 



G(E 1 ,E 2 ) = J dhdhe-^^+^Gituh). 

Equation (|E14|) then gives 

I(t) = 2eRe ]T J dEJ^S^+^g 1 ^,^ E 2 )G% R {-E 2 , E 3 ) 

ijaa' j3' 

where dE = dExdE 2 dE 3 / (2n) 5 . Using the relations [H| 

G-- = Gll — Glr, G ++ = Gll + Grl, Gll + Grl — Grr + Glr 
equation (|E16j) can be transformed to 

i+E 3 )t 



(E13) 



(E14) 



(E15) 



(E16) 



(E17) 



I(t) = 2eReJ2 J dEJ ai J* all ,e 1 ^- 



1^(E 1 ,E 2 )G^ R (~E 2 ,E 3 ) + gl R (E 1 ,E 2 )G l l + (-E 2 ,E 3 



(E18) 



We next compute the steady-state current (t — ► oo). Integrating both sides over time, we obtain the (time- 
independent) current 



/ = 2eRe ]T J ai J* a , v 



ijaa 



dE x dE 2 dE 3 



glUEuE^G^i-E.-E,) 



S ffi t (E l ,E 2 )G , L.{-E a ,-E 1 ) 



(E19) 



At steady state, since the one-particle Green's functions only depend on the difference of their arguments, we can 
make the change of variables, E\ — E 2 = E. Equation (|E19[) then gives 



/ = 2eRe £ J ai J* H , J g (gt^{E)GY R {E) + gT R (E)G\\{E) 



(E20) 



Using the relations, G*£_ — G^ + and G l [ R = —G^l Rl we can write Eq. (|E20|) as 



ii' aa' 



I = eJ2 JiaJl'f I ^ iGf^E^ilmg^m-g^iEWilmGzriE)} 



2n 



From Eqs. (|E17|) , we have -2iIm{G } = G L r + G RL - Substituting this in Eq. (|E21|) . we get 



I = eJ2 Ji«J*an> I 9l a R{E)G% L {E)-g R a L {E)Gi R {E) 



(E21) 



(E22) 



Since the leads are assumed to be in equilibrium, we can express their Green's functions using the fluctuation- 
dissipation relations 



g a L a R {E) = 2irif A (E)a aa {E)6 aa , 
g a R i{E) = 2iri(l- f A (E))a aa (E)5 a 



(E23) 
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where /a(E) — [1 + exp{f3(E — ^a)Y\ 1 an d a aa'(E) = — (l/7r)Im{G__(-E')} are the Fermi function and the density 
of states (DOS) for the left lead with chemical potential [ia- The current then becomes 



I = ieJ2rtv I ^ fA(E)G% L (E)-(l-f A (E))Gf R (E) 



2tt 



(E24) 



where T^, = ^2 aa , 2-KJi a J*, a ,o aa ' and the DOS of the leads is assumed to be energy independent in the relevant energy 
range (chemical potential difference of the two leads). This can also be expressed in terms of the self-energies due to 
the interaction with the left lead S Q/ g 



I = e J2 I T- ^R.{E)GiAE)~^L{E)Gl R {E) 



2tt 



(E25) 



where E%' R {E) = T^,f A (E) and T,% L (E) = r$,[l - f A (E)]. Equation (|E25| is same as Eq. {ITS]) . 

A similar expression can be obtained for the current from lead B to the system by replacing a with b in Eq. (|E24|) . 



i' = ie^r* / — f B {E)G% L {E) (1 - f B {E))Gf R {E) 



2?r 



(E26) 



At steady state / = —I 1 and the steady state current / s can be written as, 7 S = xl — (1 — x)I' for any value of x. 
Using Eqs. (|E24| and (|E26|) . we obtain 

— {[ x f A (E)Tt l ,-{l-x)f B (E)Tf l ,] 

If 

x {G% L {E) + G% R (E)) - [xTi, (1 - x)r£,}Gf R (E)} 



(E27) 



This can be simplified further by using, — 2iIm{G__} = Glr + G R l, to get 

Is = ~ieJ2 J ^{2i[xfA(E)r£,-(l-x)f B (E)r%,] 



x ImGTfiE) + [xT$, - (1 - x)r*]G%' R (E)} 



(E28) 



This is the general exact formal expression for the current that includes both the left and right lead properties. As 
a check let us assume no external bias and set, Ja = fs = /• Equation (|E28[) then reduces to 



-id 



E 



dE 



[xT£, - (1 - z)r* ][2i/(£)ImG»_(£) + Gf R {E)] 



(E29) 



At equilibrium, Glr(E) = —2if(E)ImG--(E) (which is the FD theorem) [fl , and I s vanishes, as it should. 

Assuming that the coupling elements are real, Yij = Fji, Eq. (|E26[) can be expressed in a simple form in terms of 
retarded and correlation Greens functions as 



I = ie£r£ J g [»(i - 2f A {E))\xaG^_{E) + \g-+{e) 



(E30) 



When the left and the right lead-system couplings satisfy T^, — XT^,, the current I s can be expressed solely in 
terms of the imaginary part of the retarded Greens function. Choosing x = 1/(1 + A) in Eq. (1E28|) gives 



ii' ii' 
" ii' ' ii' 



I"' / ^(f A (E)-f B (E))G*l_(E). 



This can be recast as, 



Is = 2eJ —(f B (E)-f A (E))J2\SiAE)f 



(E31) 



(E32) 
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where 

\S«>(E)\ 2 = ^« B lmG«_(E) (E33) 

ii* * ii' 

Equation (|E32I) has the LB form given in Eq. p. From Eq. (fE3"Tj) it is easy to see that the steady state current 
vanishes when the external bias is turned off (f A = /b) or one of the leads is disconnected (so that T A — or T B = 0). 
We note that the condition T A , = XT 1 ?, trivially holds when the leads are coupled to a single system orbital. 

When the system is modeled as non-interacting electron system and the system-lead interactions are small, we can 
use 

G i UE)=5 u , 1 ^ T - (E34) 
where a is the single electron energy. Substituting this in Eq. (|E31|) then gives 

T ° = 2e Ef5fls " /Afe)). (E35) 

l ii ' ii 

This coincide with the quantum master equation appraoch[86j. 

In general, the steady state current is given by Eq. (|E25|) which involves the Green's functions Glr and Grl- 
These can be calculated by solving the matrix Dyson equation self-consistently. When the coupling to the left and 
the right leads are proportional, T A — XT B , Eq. (|E25[) can be recast in the simpler form, Eq. (|E31[) . which involves 
properties of both the leads. Finally, when the molecule is modeled as noninteracting electron system, the current is 
simply the sum of independent contributions from each orbital (Eq. (|E35|) ). 

APPENDIX F: CONNECTING EQ. (|T49|) WITH EQ. (fT50|) 

The time evolution operator Q 1 in Eq. (|149[) is given bv[8l| 

G I {X,A,t) = fexp^-iV2 J dryi(7a(i);r)| (Fl) 



where 



vL( 7a (i);r) = -L[ik( 7£ ( T ) ;r )_ V R ( lR (r);r)] (F2) 



VMt);t) = e-^H' ABa (t) + h.c. (F3) 

Note that for ~/ a (t) = 0, ^j(7a(t);t, to) is simply the time evolution operator for the density matrix. The time 
dependence in Eq. (|F1[) is due to the interaction picture with respect to Hamiltonian Ho — Ha + Hb and H' ABa — 
J2ab J ab(^l^b) a - The parameters j a (t) = 9(t)j a with 

1L=T + 1 7a = r _2. (F4) 

-f a is an auxiliary field in Liouville space which modifies the tunneling between the two sub-systems and acts differently 
on the ket and the bra of the density matrix. A similar approach was used [2l[ to extend the well known GW equations 
22, 

[H, [3 to non-equilibrium systems. The GF, Z(X,A,t) can be computed perturbatively in the tunneling matrix 
elements, Jab- Since (V—(t)) = (particle conservation in individual noninteracting systems), to lowest (second) 
order we obtain, 

Z(\,A,t) = - [ f dTidr2(fVL(7 Q (ri);ri)VL(7 Q (r 2 );r 2 )) (F5) 
J to Jt 

By substituting Eqs. (|F2[) and (|F3p for V_ (r) and noting that (H ABa{ T i)H ABpfa)) vanishes (particle conservation; 
trace is over the product states), we obtain Eq. (|150p . 
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APPENDIX G: SELF-ENERGY FOR ELECTRON-ELECTRON INTERACTIONS 



The e — e interactions in the system are represented by the Hamiltonian. 

Hint = ^ Vijkl$i$j4>k'4>l 
ijkl 



(Gl) 



with 



Vijki = e 



dxdx'ip* (x)tfk (x) -. -rip* (x')ipi (x') 

\x — x I J 



(G2) 



and ipi is the basis set wavefunction at site i. Note that Viju — Vjuk- In this section all the superoperator indices 
shall correspond to 4- and — operators. 

The corresponding superoperator in the PTSA representation is given by 



Lint — Vijki 
ijkl af3a f 0' 



P+l, 



r,t „r,t 



(G3) 



where p is the number of "minus" operators in the second term inside the bracket. Substituting this in Eq. ()138|1 and 
expanding the exponential, the self-energy Yi a p can be evaluated perturbatively in e — e interaction, Vijki ■ To first 
order we get, 



(G4) 



ct' (3' i'k' 



where 



^S4 = (^-^)[i+(-ir +l ] 



(G5) 



with p the number of "minus" indices k', K, rj'-q of W . The electron-phonon self-energy can be calculated similarly |2£ 



APPENDIX H: SELF-ENERGY FOR A MOLECULAR- WIRE CONNECTED TO FREE ELECTRON 

LEADS 

Here we compute the Green's functions for a molecule attached to two leads. The leads are assumed to be non- 
interacting electron system which remain at equilibrium with their respective chemical potentials. We also ignore 
the many-body interactions in the molecule. For this model, the interaction with the leads can be computed exactly 
within the Green's function approach. The total Hamiltonian is given by, Eq. (|152| . 

The EOM for the Green's function G^(t, t') is obtained as in Eq. (|E5)l . 



Using the EOM for superoperator ipi a in the Heisenberg picture 

. d 



(HI) 



l-Q^Pia{t) = eii> la {t) + ^ Jix^xa, (H2) 

Eq. (|Hip can be expressed as (repeated indices are summed over), 

4 " £l ) G ^ (< ' = 5{t - tl)5ap5ij + W^M*'*' 0- (H3) 
In Eq. (|H3p . the self-energy (a) due to molecule- lead couplings is defined as 

^,{U' l )G k ^{t' l ,t') = -i Ji*(f#*a(t)ft jfi Vj). (H4) 

x£A,B 
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Note that in Eq. (|H4|) on the rhs the self-energy is given in terms of the cross-correlation Green's function defined 
in the joint space of the lead and the molecule. We wish to express it in using the Green's functions of the individual 
systems non-perturbatively. For this purpose, starting from the EOM for the lead superoperator tp xa , we construct 
the EOM for the joint Green's function. 

(i^ - <*) G%(t,t') = E J^G%(t,t'), (H5) 

^ * k 

This can also be recast as 

G%{t,t') = Y J 9 x a %{t,t'i)J*>kG k ^{t' u t') (H6) 

kx' 

where g^p is the Green's function for the free (in absence of the molecule) leads defined as, 

Substituting Eq. (|H6|) in (|H4[) . and using the identity G l ^,(t, t' x)G^, p(t' \, t') = 5 a p5ik , the self-energy can be written 

as 

Egj&O = -iJ2 J ^(t,t')J^- (H8) 

x,x' 

Since is the Green's function for the independent leads, it is diagonal in x,x', g^ = 8 X x'gap- Transforming to 
frequency domain, Eq. (|H8p then gives, 

Kp( E ) = -* E Ji*9l%E)J xj . (H9) 

xeA,B 

The self-energy can be finally obtained using the expressions for the lead Green's functions gLR. and gu,L given in 

Eqs. (|E23|) . The Green's function g^L = g-- + gLB. [Eq. (|E17|) ]. where g can be obtained from Eq. (|101[) . and 

9rr — 9ll + 9rl — 9lr [Eq. (|E17|) ]. This defines all elements of the self-energy matrix which appears in Eqs. (|119[) 
and (|E55| . 
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